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We explore the distribution of paths followed in fluctuation-induced switching between coexisting 
stable states. We introduce a quantitative characteristic of the path distribution in phase space that 
does not require a priori knowledge of system dynamics. The theory of the distribution is developed 
and its direct measurement is performed in a micromechanical oscillator driven into parametric 
resonance. The experimental and theoretical results on the shape and position of the distribution 
are in excellent agreement, with no adjustable parameters. In addition, the experiment provides the 
first demonstration of the lack of time-reversal symmetry in switching of systems far from thermal 
equilibrium. The results open the possibility of efficient control of the switching probability based 
on the measured narrow path distribution. 
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I. INTRODUCTION 

Fluctuation phenomena in systems with multiple sta- 
ble states have long been a topic of intense research in- 
terest. When the fluctuation intensity is small, for most 
of the time the system fluctuates about one of the stable 
states. Switching between the states require large fluc- 
tuations that would allow the system to overcome the 
activation barrier. Such large fluctuations are rare. How- 
ever, they lead to large changes in the system behavior. 
Fluctuational switching between coexisting states plays 
a crucial role in a variety of systems and phenomena in- 
cluding nanomagnets [l|, Josephson junctions 0], protein 
folding Q, and chemical reactions. 

A detailed theory of switching rates was first developed 
by Kramers [H . The analysis referred to systems close to 
thermal equilibrium, where the switching rate is deter- 
mined by the free energy barrier between the states. The 
concept of this barrier is well understood and the bar- 
rier height has been found for many models. However, in 
recent years much attention has been given to switching 
in systems far from thermal equilibrium. Examples in- 
clude electrons Q and atoms [a 01 in modulated traps, 
rf-driven Josephson junctions J8L [9| and nano- and mi- 
cromechanical resonators flfl fill |T2|. Apart from fun- 
damental interest, switching in modulated systems is im- 
portant for many applications, quantum measurements 
being an example [g, d Il3l . [l4| . Nonequilibrium systems 
generally lack detailed balance, and the switching rates 
may not be found by a simple extension of the Kramers 
approach. 

A basic, although somewhat counterintuitive, phys- 
ical feature of large infrequent fluctuations leading to 
switching is that, in such fluctuations, the system is most 
likely to move along a certain path in its phase space. 
This path is known as the most probable switching path 
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(MPSP). For low fluctuation intensity, the MPSP is ob- 
tained by the solution of a variational problem. This 
problems also determines the switc hing activation bar- 
rier [H, M, 03, [H, EI HI HI, HI, H, H| . Despite its 
central role for the understanding of fluctuation-induced 
switching and switching rates, the idea of the MPSP has 
not been tested experimentally in multivariable systems. 
The question of how the paths followed in switching are 
distributed in phase space has not been asked either. 

Perhaps closest to addressing the above issues was the 
experiment on dropout events in a semiconductor laser 
with optical feedback [2(| . In this experiment the switch- 
ing path distribution in space and time was measured 
and calculated. However, the system was characterized 
by only one dynamic variable, and thus all paths lie on 
one line in phase space. In another effort, electronic cir- 
cuit simulations [271 ] compare the distribution of fluctu- 
ational paths to and relaxational paths from a certain 
point within one basin of attraction; the data refer to the 
situation where switching does not occur. The methods 
[H H3 do not apply to the switching path distribution 
in multivariable systems, as explained below. 

In the present paper we introduce the concept of 
switching path distribution in phase space and the quan- 
tity that describes this distribution, calculate this quan- 
tity, and report the direct observation of the tube of 
switching paths. A brief account of the results was pub- 
lished in Ref. [28|. The experimental and theoretical re- 
sults on the shape and position of the path distribution 
are in excellent agreement, with no adjustable parame- 
ters. The results open the possibility of efficient control 
of the switching probability based on the measured nar- 
row path distribution. 

In Sec. [II] we provide the qualitative picture of switch- 
ing and give a preview of the central theoretical and 
experimental results. Sec. IIIII presents a theory of the 
switching probability distribution in the basins of attrac- 
tion to the initially occupied and initially empty stable 
states as well as some simple results for systems with de- 
tailed balance. In Sec.[lV]the system used in the experi- 
ment, a micromechanical torsional oscillator, is described 
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and quantitatively characterized. Section [V] presents the 
results of the experimental studies of the switching path 
distribution for the micromcchanical oscillator, with the 
coexisting stable states being the states of parametri- 
cally excited nonlinear vibrations. Generic features of the 
distribution are discussed, and the lack of time-reversal 
symmetry in switching of systems far from thermal equi- 
librium is demonstrated for the first time. Section IVII 
contains concluding remarks. 

II. QUALITATIVE PICTURE AND PREVIEW 
OF THE RESULTS 

We consider a bistable system with several dynamical 
variables q = (qi, q^). The stable states A\ and A 2 
are located at (\a 1 and q^ 2 respectively. A sketch of the 
phase portrait for the case of two variables is shown in 
Fig. [TJ For low fluctuation intensity, the physical picture 
of switching is as follows. The system prepared initially 
in the basin of attraction of state A\, for example, will 
approach over the characteristic relaxation time t r 
and will then fluctuate about q^ . We assume the fluctu- 
ation intensity to be small. This means that the typical 
amplitude of fluctuations about the attractor (the char- 
acteristic diffusion length) Ijj is small compared to the 
minimal distance Ag between the attractors and from 
the attractors to the saddle point qs. 

Even though fluctuations are small on average, occa- 
sionally there occur large fluctuations, including those 
leading to switching between the states. The switching 
rate W12 from state A\ to A2 is much less than the recip- 
rocal relaxation time i~ , that is, the system fluctuates 
about A\ for a long time, on the scale of t r , before a 
transition to Ai occurs. In the transition the system 
most likely moves first from the vicinity of q J 4 1 to the 
vicinity of qs- Its trajectory is expected to be close to 
the one for which the probability of the appropriate large 
rare fluctuation is maximal. The corresponding trajec- 
tory is illustrated in Fig. [TJ From the vicinity of qs 
the system moves to state A2 close to the deterministic 
fluctuation-free trajectory. These two trajectories com- 
prise the MPSP. 

For brevity, we call the sections of the MPSP from q^ x 
to qs and from qs to (\a 2 t ne downhill and uphill tra- 
jectories, respectively. The terms would literally apply 
to a Brownian particle in a potential well, with A\$, cor- 
responding to the minima of the potential and S to the 
barrier top. 

We characterize the switching path distribution by the 
probability density for the system to pass through a point 
q on its way from A\ to A2, 

PM<l,t)= dq f p(qf,tf,(i,t\q ,to). (1) 

Here, the integrand is the three-time conditional proba- 
bility density for the system to be at points q^ and q at 
times tf and t, respectively, given that it was at qo at 
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FIG. 1: (Color online) Phase portrait of a two- variable system 
with two stable states Ai and Ai. The saddle point S lies on 
the separatrix that separates the corresponding basins of at- 
traction. The thin solid lines show the downhill deterministic 
trajectories from the saddle to the attractors. A portion of 
the separatrix near the saddle point is shown as the dashed 
line. The thick solid line shows the most probable trajectory 
that the system follows in a fluctuation from A\ to the sad- 
dle. The MPSP from Ai to A2 is comprised by this uphill 
trajectory and the downhill trajectory from <S to A2- The 
plot refers to the system studied experimentally, see Sec. |V] 

time to- The point qo lies within distance ~ of q^ 
and is otherwise arbitrary. Integration with respect to q^ 
goes over the range of small fluctuations about q^ 2 ; 
the typical linear size of this range is Ip. 

We call P12 (q, t) the switching probability distribution. 
Of utmost interest is to study this distribution in the time 
range 

W21 > tf - t, t - t > t r . (2) 

Here, t r is the Suzuki time [29(. It differs from t r by a 
logarithmic factor ~ log[Ag/Zo]. This factor arises be- 
cause of the motion slowing down near the saddle point. 
The time t r is much smaller than the reciprocal switch- 
ing rates, and the smaller the fluctuation intensity the 
stronger the difference, because the dependence of Wij on 
the fluctuation intensity is of the activation type. If the 
noise causing fluctuations has a finite correlation time, t r 
is the maximum of the Suzuki time and the noise corre- 
lation time. 

For t — to 3> t r , by time t the system has already 
"forgotten" the initial position qo. Therefore the distri- 
bution p(qf,tf] q, i|qo, to), and thus P12, are independent 
of qo,^o- On the other hand, if the system is on its way 
from A\ to A2 and at time t is in a state q far from the 
attractors, it will most likely reach the vicinity of A2 over 
time < t r and will then fluctuate about q^ 2 • This will 
happen well before the time tf at which the system is 
observed near A 2l and therefore p\2 is independent of tf. 

It is clear from the above arguments that, in the time 
range ((U), the distribution pi2(q, t) for q far from the 
attractors is formed by switching trajectories emanating 
from the vicinity of A\ . It is gives the probability density 
for these trajectories to pass through a given point q 
at time t. In other terms, the distribution pi2(q, t) is 
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formed by the probability current from Ai to A2 and is 
determined by the current density. 



The shape of the switching probability 
distribution 



We show in Sec. Hill that j>i 2 (q,t) peaks on the MPSP. 
The peak is Gaussian transverse to the MPSP for |q — 

q^i, 2 l, |q- qsl > Id, 



Pi 2 (q,*) = Wutr^fnJZ^exp --£±_Q£ 



1 



(3) 



where £|| and £j_ are coordinates along and transverse 
to the MPSP, and v{(,\\) is the velocity along the MPSP. 
The matrix elements of matrix Q = <2(£ii) are oc , and 

Z = [(2vr) JV - 1 / dot Q] 1 ' 2 . It follows from Eq. © that the 
overall probability flux along the MPSP is equal to the 
switching rate, 



d£-LPi2(q,tM£||) = W 12 . 



We have observed a narrow peak of the switching 
path distribution in experiment. The results are shown 
in Fig. [21 They were obtained using a micro-electro- 
mechanical torsional oscillator described in Sec. HVl The 
path distribution displays a sharp ridge. We demonstrate 
that the cross-section of the ridge has Gaussian shape. As 
seen from Fig. [21 the maximum of the ridge lies on the 
MPSP which was calculated for the studied system. 
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FIG. 2: (Color) (a) Switching probability distribution in a 
parametrically driven microelectromechanical oscillator. The 
probability distribution p\2{X,Y) is measured for switching 
out of state Ai into state Ai. (b) The peak locations of the 
distribution are plotted as black circles and the theoretical 
most probable switching path is indicated by the red line. All 
trajectories originate from within the green circle in the vicin- 
ity of A\ and later arrive at the green circle around Ai. The 
portion of the distribution outside the blue lines is omitted. 

Equation ([3]) is written for a generally nonequilibrium 
system, but the system is assumed to be stationary. In 
the neglect of fluctuations its motion is described by 
equations with time-independent coefficients. In this case 



Pi2(q, t) is independent of time t. A different situation 
may occur in periodically modulated systems. In such 
systems, attractors are periodic functions of time. If the 
typical relaxation time is smaller than or of the order of 
the modulation period, generally there is one MPSP per 
period. Then pi2(q, t) is also a periodic function of time. 
We will not consider this case in the present paper. 



B. Comparison with the prehistory distribution 

The distribution of fluctuational paths was studied ear- 
lier in the context of the "prehistory problem" [3(| ■ In 
this formulation one is interested in the paths to a certain 
state q/ that is far from the initially occupied attractor. 
The distribution of these paths ph is given by the prob- 
ability density to have passed a point q at time t given 
that the system is found at q/ at a later time t / whereas 
initially at time to it was at point qo near attractor A\, 



Ph(<ht\qf,t f ) 



p(q/,*/;q,*|qo,*o) 
p(q/»*/;qo)*o) 



(4) 



The prehistory distribution ([4]) and its generalizations 
were analyzed in a number of papers H,[23,[3l|,[l,[3l. 
However, the problem of paths that lead to switching 
between the states was addressed only for a stationary 
system with one dynamical variable [261 ]. In this case, 
the system must pass through all the intermediate points 
between the two states during a switch. For systems 
with more than one dynamical variable, the formula- 
tion [2^, [27], HH no longer applies because it cannot be 
known in advance through what points the system will 
pass in switching. The aforementioned formulation does 
not work even for one-variable periodically modulated 
systems, since the distribution ([4]) depends not only on 
the position of point qr, but also on the time tf when 
this point is reached [34| . 

In contrast, the distribution pi2(q, t) is defined in such 
a way that it is independent of the final point qy and of 
the time tf of reaching it . The definition does not impose 
any constraint on paths except that they lead to switch- 
ing between the attractors. Therefore the introduction of 
the function pi2(q, t) is essential in studying the switch- 
ing path distribution for multivariablc systems. 



III. THEORY OF THE SWITCHING PATH 
DISTRIBUTION 

A. The model of a fluctuating system 

We derive Eq. (2) for a system described by the 
Langevin equation of motion 

q = K(q) + f(t), (f n (t)f m (t')) = 2DS nm S(t - t'). (5) 

Here, the vector K determines the dynamics in the ab- 
sence of noise; K = at the stable state positions q^ , 
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(\A 2 and at the saddle point q,g. We assume that c\$ lies 
on a smooth hypersurface that separates the basins of 
attraction of states A\ and A 2 , cf. Fig. [TJ The function 
f(t) in Eq. © is white Gaussian noise; the results can 
be also extended to colored noise. The noise intensity D 
is assumed small. The dependence of the switching rates 
Wnm on D is given by the activation law, log W nm oc D^ 1 
M EE E3, El EE HE E3 . This is also the case for noise- 
driven continuous systems, cf. Ref. HE HE H3 and papers 
cited therein. There exists extensive literature on nu- 
merical calculations of the switching rate and switching 
paths, cf. Ref. [38, 39, 4(j l4ll l42l and papers cited therein. 

In the model ([5]), the characteristic relaxation time t r 
and the characteristic diffusion length Id are 



t r = max I Re Xk\ 

k 



Id = {Dt r 



:V2 



(G) 



where Aj, are the eigenvalues of the matrix dK m /dq n cal- 
culated at qAj , Ha 2 and qg. 

For a white-noise driven system the three-time 
probability distribution /o(q/, tf, q, t|qo, to) in Eq. ((T|) can 
be written as a product of two-time transition probability 
densities, 



p(q/,*/;q>*|qo,*o) = p(q/,*/|q,*)p(q>*|qo,*o 



(7) 



which simplifies further analysis. The analysis is done 
separately for the case where the observation point q 
lies within the attraction basins of the initially empty 
attractor A 2 and the initially occupied attractor A\ . 



B. Switching probability distribution in the 
initially unoccupied basin of attraction 

We start with the case where the observation point q 
lies in the basin of attraction of the initially empty state 
A2 far from the stationary states, |q — qs|, |q — q^i 2 1 ^ 
Id ■ For weak noise intensity, the system found at such q 
will most likely approach q^ 2 over time t x moving close 
to the noise-free trajectory q = K and will then fluctuate 
about qA 2 ■ Therefore, for q/ not far from the attractor 
A 2 , i.e., |q/-<3U 2 | < Id, we have p(qj,t f \q,t) w /0 2 (q/)- 
Here, p 2 {<\f) is the stationary probability distribution 
in the attraction basin of A 2 in the neglect of A 2 — > 
A\ switching. In its central part it has the form of a 
normalized Gaussian peak centered at q^ , with typical 
size Id- Then, from Eq. |T]) 

Pi2(q,i) = p(q,t|qo,*o)- 

The analysis of the transition probability density 
p(q, t|qo, to) in this expression is simplified by two obser- 
vations. First, for time t in the range W^ 2 3> t— to 3> t r , 
there is a probability current from attractor A\ to A 2 . 
This current gives the switching rate W\ 2 , as found by 
Kramers The current density far from A 2l i.e., for 
q~ q^ 2 1 ^ Id, is independent of time and is determined 
by the stationary Fokkcr-Planck equation 



The second observation is that, for both white and 
colored Gaussian noise, in switching the system is most 
likely to go close to the saddle point [IE [43|. Having 
passed through the region near the saddle point the sys- 
tem moves close to the deterministic downhill trajectory 
from S to A 2 , cf. Fig. [TJ This trajectory is described by 
equation q = K and gives the MPSP in the basin of at- 
traction of A 2 . We are interested in finding p(q, t|q , to) 
for q close to this trajectory. The broadening of the dis- 
tribution is due to diffusion, which shou ld g enerally make 
it Gaussian in the transverse direction [44( . 

We parameterize the deterministic section of the 
MPSP by its length £|| counted off from qs and intro- 
duce a unit vector £\\ along the vector K on the MPSP 
and N — 1 vectors £j_ perpendicular to it. The velocity 
on the MPSP is v = w(^) = =0). Of interest 

for our analysis are the values of of the order of the 
width of the path distribution transverse to the MPSP, 
which is given by the diffusion length, i.e., < Id- 
We assume to be small compared to the radius of 
curvature \d^/d^ 

Equation {8} can be solved near the MPSP by changing 
to variables £|| , expanding K to first order in and 
replacing — > d^ ± . One then obtains for p(q, t\q ,to) = 
pi2(q, t) expression l[8]). with matrix Q given by equation 



v-j— + k'Q 



Qk + 2Q 2 D = 0. 



(9) 



Here, = dK^/d^^v, with the derivatives evaluated 
for £± = 0; the subscripts fj,, v = 1, . . . , N — 1 enumerate 
the components of £j_ and the transverse components of 
K in the co-moving frame. Equation @ can be reduced 
to a linear equation for Q -1 . From Eq. ([9]), the matrix 
elements Q^^ ocl/D. Therefore the width of the switch- 
ing probability distribution (2) is oc Id, as expected from 
qualitative arguments. 



C. Switching probability distribution in the 
initially occupied basin of attraction 

The case where the observation point q lies in the basin 
of attraction of the initially occupied state A\ is some- 
what more complicated. Here, too, the two probability 
densities in the right-hand side of Eq. {ZJ are independent 
of time t for |q — q^i 2 |j |q _ qsl ^ Id- But in contrast to 
the previously studied region, none of them is known in 
advance. They have to be found from the Fokkcr-Planck 
equation (|5| for p(q, t|qo,io) and the backward equation 
for p(q/,t/|q, t), 

(Ka q + D^)p(q / ,t / |q,t) = 0. (10) 

We seek the solutions of Eqs. (jSJ and (TlO|) in the cikonal 
form, 



D$?]p(q,t|cto,to)=0. 



(8) 



/9(q,t|qo,*o) = exp[-S F (q)/D], 
p(q/,*/|q,*) = cxp[S , B (q)/-D]p2(q/)- 



111) 
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The functions 5f and 5b can be written as power series 



in the noise intensity D, with 5f, b 
To the lowest order in D we have 
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(o) 

F, B 



-DS, 



(i) 

F, B 



# (q, d q 5 F 0) 



0, 



i?(q,p) = p 2 + pK(q). (12) 



Equation (fT2"|) has the form of a Hamilton- Jacobi equa- 
tion for an auxiliary particle with coordinate q and mo- 
mentum p. This particle moves with energy H = 0. 
The functions 5 F B (q) are mechanical actions. Subscript 
F refers to motion of the auxiliary particle to point q 
from the vicinity of Ai, as it is clear from Eq. (fTTj) . 
From Eq. ([2|) , this motion takes time that largely exceeds 
t r . Using condition H = 0, one can therefore associate 

5 F (q) with the mechanical action for reaching q from 
flAx', the motion formally starts at t — > -co from q^i, 
with momentum p = [151 ]. 

Subscript B in Eq. (TT2]) refers to the auxiliary Hamil- 
tonian particle that moves from q further away from at- 
tractor A\ . In this motion the original system goes close 
to the saddle point, and so should the auxiliary particle, 
too. The perturbation theory that underlies Eq. (|12[) ap- 
plies where the particle is approaching the saddle point, 
but has not gone beyond it. Indeed, for H = the par- 
ticle approaches the saddle point asymptotically, for in- 
finite time. Therefore 5p°^ is the mechanical action for 
reaching q^ from q. 

From Eqs. 0, fir]), the MPSP inside the basin of at- 
traction of the initially occupied state corresponds to the 

maximum of (q) — 5 F °' ) (q) and thus is determined by 
equation 



d q 5 F 



(o) 



a q 5 



(0) 



(13) 



The MPSP is thus given by the hctcroclinic Hamilto- 
nian trajectory that goes from the state (q^ , p = 0) to 
(qs,P = 0). 

To find 5 F ° B (q) close to the MPSP it is convenient 
to switch to a co-moving frame on the MPSP 
From Hamiltonian (|12p . the longitudinal direction £\\ and 
the velocity on the MPSP are given by expression 



2d q 5 F 0) (q) + K(q)=t,(£||)£||, 



(14) 



where the left-hand side is calculated for £± = 0. 
[Eq. (TT4"]) applies also if we use 5 B 0) instead of 5 F 0) ]. Note 

that the MPSP direction £\\ is not along the velocity of 
the original system in the absence of noise K, in the gen- 
eral case of a system lacking detailed balance. 

Close to the MPSP we can expand 5 F ° B and K 
in £ ± . From Eqs. Q, HJ) and (T3J), Pi2(q,*) oc 
cxp(— £j_Q£_l/2), as in Eq. ([3]). The matrix Q is ex- 
pressed in terms of the actions 5 F ° B close to the MPSP 
as 

Q = Qf — Qb, 

(Qf,b)„„ = D- 1 d 2 S ( p %/d^d^ v , (15) 



with the derivatives calculated for £j_ = 0. From the con- 
dition that pi2(q, t) be maximal on the MPSP it follows 
that matrix Q is positive definite. 



1. The prefactor 

Interestingly, the prefactor in p\2 (q, t) can be ex- 
pressed explicitly in terms of the velocity i>(£||) and the 
matrix Q. Formally, the prefactor is determined by the 
terms 5 F B in Eq. pip . The equations for them follow 



from Eqs. © and (T 



23 q 5 F 



(o) 



K 



<9 q 5p 



a q K 



9 q 5 F 



(0) 



0. 



2d q 5 B 



(o) 



K 



<9 q 5 B 



(i) 



(16) 



From Eqs. (Tj"4")) . (fTB]), to leading order in £j_ we have 



vd e , 



Tr 



D 



(Qf + Qb) 



0, 



5« = 5 F 1) -5 B 1) , 



(17) 



where, as before, = dK^/d^u with the derivative 
calculated for £± = 0. 

On the other hand, by expanding in Hamilton- Jacobi 
equations (|12p for 5 F ° B near the MPSP to second order 
in £± and taking into account the relation between the 

derivatives of 5 F 0) and 5^ 0) on the MPSP 0, we 
obtain an important relation 



vd^Q + 2D 
From this equation 
Tr 



Qf 



Qk = 0. 



k + £>(Qf + Qb) = --«0 £|| Trlog< 



By substituting this relation into Eq. I|17p we obtain 

S<» (C|| ) = log «($| ) - ^Tr log Q(C|| ) + log & , (18) 

where we explicitly indicate that v, and Q are func- 
tions of the distance, along the MPSP; C\ is a con- 
stant of integration. 

Equations |T]), (fTTj) , (fl5|) . and (JTHJ) lead to expression 
Q for the switching probability distribution. Note that, 
from Eq. (fTS"|) . inside the initially occupied basin of at- 
traction the width of the peak of the distribution trans- 
verse to the MPSP is ~ Irj oc D 1 / 2 . The distribution 
describes a stationary probability current. This current 
is the same in the basins of attraction of states A\ and 
A<2,. In obtaining Eq. © from Eq. (fT5|) we found C\ from 
the condition / d£_i_Pi2 (q, t) = W\2- 

From conservation of the stationary probability current 
it follows that the distribution p\2 (q, t) should sharply 
increase near the saddle point. Indeed, the velocity 
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= for q = q$. The current close to qs is due 
to diffusion. In the general case of nonequilibrium sys- 
tems the shape of the switching probability distribution 
near the saddle point is complicated; its analysis is be- 
yond the scope of this paper. 



D. Switching probability distribution for systems 
with detailed balance 



An explicit solution for pi2(q, t) near the saddle point 
can be obtained for systems with a gradient force K = 
— dqU(q). Such systems have detailed balance. The up- 
hill section of the MPSP is literally the uphill path that 
goes from the local minimum of the potential t/(q) at Ai 
to the saddle S and is given by equation q = <9 q [7 (q) [45[ 
[this can be seen from Eq. (1121) 1. In contrast to systems 
without detailed balance [43|, [46| , for smooth U (q) the 
MPSP near the saddle point is described by an analytic 
function of coordinates and £m is perpendicular to the 
separating hypersurface. 

The quasistationary solution of the forward Fokkcr- 
Planck equation (jSJ) near a saddle point has been known 
since the work of Kramers Q and Landauer and Swanson 
[47l ] . The backward equation (fpQ|) can be solved similarly 
by expanding K to first order in q — qs and by using the 
condition that deep inside the basin of attraction of the 
initially empty state A2 we have p(q/,t/|q, t) « /^(q/)- 
The solution has the form 

p(q/,*/|q,*) ~ 5P2(q/) l + erf(|||) 

fll^AiiAD) 1 /^,,^). (19) 

Here, erf (a;) is the error function, £115 is the position of 
the saddle point on the MPSP, and A|| is the curvature 
of the potential U(q) at the saddle point in the steepest 
descent direction £|| , C/(q) ~ — A|| (£|| — £|| s) 2 /2 for £j_ = 
and small |£|| -£||s|. 

Equation (flT))) combined with the results [J, |47| give 
expression ([3J forpi2(q, t) near qs provided one replaces 
in this expression 



— (7r/8A|,Z?) 1 / a exp(^) 



1- erf 2 (|,| 



(20) 



Equation ([20)) goes over into f _1 (£||) for £|| — £||s| ^ 
In the opposite limit, that is very close to the 
saddle point, it shows that v~ x is replaced by a factor 
(7T/8AH-D) 1 / 2 ~ t r /ln- This demonstrates that the dis- 
tribution pi2(q, t) does not diverge at the saddle point, 
but it contains a large factor D~ x l 2 . 



IV. MICROMECHANICAL TORSIONAL 
OSCILLATOR 

A. Device characteristics 

We measure the switching probability distribution us- 
ing a high-Q micro-electromechanical torsional oscillator 



(Q = 9966) driven into parametric resonance. The oscil- 
lator is shown in Fig. [31 It consists of a movable, highly- 
doped polysilicon plate (200 /mi x 200 (im x 3.5 /jm) 
suspended by two torsional rods (4 /im x 2 /im x 36 /im, 
spring constant = 3.96 x 10~ 8 Nm). There are two fixed 
electrodes on the substrate, one on each side of the tor- 
sional rod. The 2 /im gap underneath the movable plate 
is created by etching away a sacrificial silicon oxide layer. 





FIG. 3: Micromechanical torsional oscillator used for study- 
ing the switching path distribution, (a) Scanning electron 
micrograph, (b) Cross-sectional schematic. The angle 9 of 
the movable plate is controlled by the voltage applied to one 
of the fixed electrodes. Oscillation of the plate is detected 
using the other electrode. 

Torsional oscillations of the movable top plate are 
excited by applying a driving voltage Vd = Vdc + 
V ac cos (ujt) + V no i se (t) to one of the lower electrodes while 
the top plate remains electrically grounded. The driving 
frequency uj — 2cu + e is close to twice the natural fre- 
quency u>o- The dc voltage Vd c (1 V) is much larger than 
the amplitude Va C (141 mV) of sinusoidal modulation and 
the random noise voltage Vn i sc . 

Because of the applied voltage, the top plate is sub- 
jected to an electrostatic torque r = (dC / d9)V d 2 / 2 in 
addition to the restoring torque of the springs: 

9 + 2T6 + uj18= I (21) 



where 9 is the rotation angle (see Fig. [3j, T is the damp- 
ing constant, and / is the moment of inertia of the 
plate. The time-independent component of the torque 
G = — C'(9 )V^ c /2 leads to a shift of the equilibrium po- 
sition of 9 to 9 . In what follows we count 9 off from 
- 

The expansion of the torque in 9 to first order in 
Vac, Koise has the form 



f39 3 



— [k + k e cos {uit)) 9 — a 
+F cos (u>t)+N(t). (22) 

Here, k - -C"(9 )Vl/2I, a = (6 )Vl/4I, and 

[3 = —C^ (9 )V^ C /12I are the linear and nonlinear coeffi- 
cients of the electrostatic torque; because the oscillation 
amplitude remains small, we disregard terms of higher 
order in 9. 

The ac voltage leads to a time-dependent additive 
torque with amplitude F = C (9o)Vd c V :xc (t)I and also 
to modulation of the spring constant oc k e , with k e = 
—C" (9o)Vd c V ac / 1 . Since the Q-factor of the oscillator 
is high, the response at ui ~ 2u> is negligible. There- 
fore the major effect of the ac voltage is the paramet- 
ric modulation k e 9 cos (cut). The renormalization of this 
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term by the nonlinear terms in Eq. (|22|) [for example, 
oc aOujQ 2 F cos cut] is small for our device and is disre- 
garded. The term N(t) — C (9o)Vd c V no i se (t) / 1 represents 
zero-mean noise in the driving torque. This noise is Gaus- 
sian and its spectrum is flat in a broad frequency range 
that goes from zero frequency far beyond ujq] therefore 
for our purpose it can be assumed white, with intensity 
D given by (N(t)N(t')) = 2D5(t - t'). 

With Eq. (|22|) , the equation of motion for the angle 9 
counted off from 9 becomes 

9 + 2T9 + [uf + k e cos (cut)] 9 + a9 2 + /36> 3 = N(t), (23) 

where to 2 = ujq +ko (in our device the difference \ujo — <^i \ 
is small, < O.OOlwo). 

Torsional oscillations of the top plate are detected ca- 
pacitively by the other electrode. This electrode is con- 
nected to a dc voltage source through a large resistor. A 
high electron mobility transistor is placed in close prox- 
imity to the device to measure the oscillating charge on 
the detection electrode induced by motion of the top 
plate. The output of the transistor is connected to a 
lock-in amplifier referenced at half the driving frequency 
uj. For the chosen time constant of 300 [is, the measure- 
ment uncertainty is ~ 80 /irad, about 0.6% of the full 
scale in Fig. [5] and much smaller than the width of the 
path distributions. The oscillation amplitudes in-phasc 
(X) and out-of-phase (Y) with the reference frequency 
were recorded every 2 ms. All measurements were per- 
formed at 77 K and < 10~ 6 torr. 



B. Transformation to slow variables and 
parametric resonance 

Since the oscillator is strongly undcrdamped 
(r/u>i ~ 10~ 4 ) and the modulation is almost reso- 
nant (\u> — 2u)i\ <C u>), we analyze the motion of the 
oscillator in the rotating frame, with slow dimensionless 
variables q± and q 2 and dimensionless time t — » k e t/2uj 
(note that even though the oscillator has one degree of 
freedom, its motion is characterized by two dynamical 
variables). In our micromechanical device the charac- 
teristic rcnormalizcd parameter of cubic nonlincarity 
7 = (3 — (10a 2 /9uj 2 ) is negative. In this case it is 
convenient to introduce the slow variables as 



m = ^7 



2k, 
3M 
ui 2 k t 



1/2 



1/2 



91 cos "TT 



q 2 sm — 



qi sm 



92 cos 



(24) 



The variables qi and q 2 are interchanged here compared 
to Ref. |H, which referred to the case 7 > 0. 

The quadratures 91 and 92 are directly proportional to 
the signal components X and Y measured with the lock- 
in amplifier, with the proportionality constant E deter- 
mined by the measuring apparatus, 



Substituting Eq. (|24|) into Eq. ([23]) and neglecting fast 
oscillating terms, we can write the equations of motion 
for q = (91,92) in the form ([5]). The function K in di- 
mensionless time is given by 



K 



-rV 



eVg. 



(26) 



uj(2oji—uj)/k e , and g = q 4 /4 — (1- 



91 = EX, 92 = EY. 



(25) 



Here £ = k e /2ujT, fx 

fx)q 2 /2 + (1 — [i)q 2 /2, where e is the permutation tensor 
(the parameter ji is defined with the opposite sign com- 
pared to Ref. |48], again because 7 has the opposite sign 
for our system) . Equation q = K gives the downhill sec- 
tion of the MPSP of the oscillator. The uphill section of 
the MPSP can be calculated by solving the Hamiltonian 
equations of motion that follow from Eq. (fTS")) . 



C. Determination of device parameters 

We first consider motion of the device in the absence 
of fluctuations. When the amplitude of the spring mod- 
ulation is sufficiently strong (£ > 1), the oscillator re- 
sponse exhibits period doubling [49(. Oscillations arc 
induced at half the modulation frequency in a range 
close to wi. Between the two bifurcation frequencies u)bi 
(w 139318.11 rad/s) and uu b2 (k, 139384.74 rad/s) there 
exists two stable states of oscillations at frequency u/2. 
They differ in phase by 7r but have identical amplitude. 
Both states are stable solutions of Eq. (|23]) . Their basins 
of attraction in the rotating frame are separated by a sep- 
aratrix that goes through the unstable stationary state, 
which in the laboratory frame has zero vibration ampli- 
tude at frequency uj/2. The phase portrait in the rotat- 
ing frame is illustrated in Fig. [TJ The driving frequency 
is chosen to be 278639.16 rad/s for measurement of the 
switching path distribution. 

We note that parametric resonance in nano- and micro- 
electro-mechanical systems has attracted considerable at- 
tention EE IS, HI, 0, H5[. Since here we are 
interested in the studies of the principal features of 
noise-induced switching, we chose the simplest nontriv- 
ial regime where the system has only two stable states, 
which occurs for lu^i < lu/2 < u)b 2 . The modulation fre- 
quency oj is chosen to be close to 2u!ti so that the motion 
in the rotating frame is underdamped, which is advan- 
tageous for studying a generic feature of fluctuations in 
systems far from thermal equilibrium, the breaking of 
time reversal symmetry. 

Calculation of the MPSP requires a number of device 
parameters including T, u>\, the parametric modulation 
amplitude k e , and the nonlinear constant K non ii, loar = 
37/8W1 [H[. These parameters are obtained from the 
linear and nonlinear responses of the device. When the 
device is resonantly driven with small amplitude at fre- 
quency close to cJi, it responds as a harmonic oscillator. 
From the resonance line shape (Fig. [4j, T and u>\ are 
determined to be 6.99 rad/s and 139352.118 rad/s re- 
spectively. The remaining two parameters are extracted 
from the parametric resonance of the oscillator for tu close 
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FIG. 4: (Color online) Harmonic and parametric resonances 
of the micromechanical torsional oscillator. For resonant driv- 
ing (solid circles), the oscillation amplitude is plotted as a 
function of the oscillation frequency. The thin line is a fit to 
the harmonic oscillator response. It gives device parameters T 
and u>i . For parametric resonance (hollow squares) , the driv- 
ing frequency is twice the oscillation frequency. The fit (thick 
line) yields Knonimear and the effective parametric modulation 
amplitude fc e . 

to 2u>i. Specifically, the parametric modulation ampli- 
tude k e is determined from the bifurcation frequencies 
Lu b i^ = 2u>i T ijJ p , where uj p = (fc 2 - (AloiT) 2 ) 1 ^ 2 /2loi. 
This gives k e = 1.94 x 10~ 7 s~ 2 . The nonlinear parame- 
ter K non i incar (1.08xl0 6 s _1 ) is obtained from the propor- 
tionality constant between the square of the parametric 
oscillation amplitude 9a and the detuning from bifurca- 
tion frequency close to the bifurcation frequency seen in 
Fig.S 

®A = { u ~ w &2)/2ftnonlmcar. (27) 

Using these measured device parameters, the dimen- 
sionlcss constants contained in K in Eq. (|26[) and in 
Eq. (Jini) can be directly calculated to be E = 176.349, C 
= 4.968, and fi =-0.9367. The theoretical optimal escape 
path in Fig. [2] is calculated with the above parameter 
values. No adjustable parameters are used. 

V. SWITCHING PATH DISTRIBUTION: 
EXPERIMENT 

A. Measured switching path distribution 

When white noise is added to the excitation voltage, 
the system can occasionally overcome the activation bar- 
rier and switch from one stable state to the other. The 
noise intensity is chosen to ensure that the mean resi- 
dence time in each state (~ 10 s) is much larger than 
the relaxation time (t r ~ 1 s) of the system. Transitions 
are identified when the oscillator begins in the vicinity of 
A\ (within the left green circle fii in Fig. [5^,) and subse- 
quently arrives at state Ai (within the right green circle 



) ■ Figure [2] shows the switching probability distribu- 
tion derived from more than 6500 transitions. While in 
each transition the system follows a different trajectory, 
the trajectories clearly lie within a narrow tube. 

The maximum of the distribution gives the MPSP. In 
Fig. f2p, the location of this maximum is plotted on top 
of the MPSP obtained from theory. The oscillator is 
underdamped not only in the laboratory frame, but also 
in the rotating frame. Therefore both the uphill and 
downhill sections of the MPSP are spirals. On the uphill 
section, the MPSP emerges clockwise from A\ and spirals 
toward the saddle point at the origin. Upon exit from 
the saddle point, it makes an angle and, on the downhill 
section, continues to spiral clockwise toward A^. 

There is excellent agreement between the measured 
peak in the probability distribution and the MPSP ob- 
tained from theory. There are no adjustable parame- 
ters since all device parameters are accurately determined 
from the harmonic and parametric resonances of the os- 
cillator without noise in the excitation as described in 
the previous section. 

Close to the stable states the peaks of the distribution 
at successive turns of the MPSP overlap, preventing the 
accurate determination of the MPSP. The plot in Figs. [2^ 
and f5p has excluded the portions of trajectories prior 
to escaping from the initial state A\ and upon arriving 
at the final state A2, which are bound by the two blue 
lines. Such cutoff also eliminates the large peaks of the 
distribution centered at Ai and A2, which arise because 
the oscillator spends most of its time fluctuating about 
A\ and A2- These peaks are not relevant to switching 
dynamics. 

Figure [5] compares the measured and predicted veloc- 
ity along the MPSP. Here, again, the good agreement 
is demonstrated with no adjustable parameters. As ex- 
pected, the measured velocity decreases near the saddle 
point, £|| = 0. However, it does not become equal to zero, 
in agreement with the argument that the total probabil- 
ity current remains constant. Motion near the saddle 
point is dominated by diffusion. 
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FIG. 5: Measured averaged velocity along the MPSP (cir- 
cles) and the velocity predicted by theory (line) . The velocity 
decreases to zero at the saddle point (fn = 0). 
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B. Generic features of the switching path 
distribution 

The switching probability distribution in our multi- 
variable system displays important generic features. Fig- 
ure [5£i shows the distribution cross-section along the pur- 
ple line transverse to the MPSP in Fig. [2}d. It is well- 
fitted by a Gaussian. Gaussian distributions with dif- 
ferent height and area are observed also in other cross- 
sections except close to the saddle point. Figure [5b plots 
the area under the Gaussian distribution versus the recip- 
rocal measured velocity on the MPSP, for different cross- 
sections. The linear dependence agrees with Eq. and 
indicates that the probability current from the initially 
occupied attractor to the empty one is constant. This 
current gives the switching rate W\2 

We find that the probability current concentrates 
within a narrow tube deep into the basins of attraction 
of A\ and A2. In the basin of attraction to A2 but not 
too close to A2, much of the probability distribution car- 
ries the switching current. However, the overall quasis- 
tationary probability distribution deep inside the basin 
of attraction of A\ is largely associated with fluctuations 
about A\ that do not lead to switching. The part of 
the distribution responsible for the switching current is 
an exponentially small fraction of the total distribution. 
Nevertheless our formulation makes it possible to single 
out and directly observe this fraction. 




e(mrad) 1/Area(rad) 

FIG. 6: (a) The cross-section along the purpleline in Fig. [2] 
transverse to the MPSP. The solid line is a Gaussian fit. (b) 
Velocity on the MPSP vs. inverse area under cross-sections 
of the switching probability distribution. The solid line is a 
linear fit forced through the origin. 



The slowing down near the saddle point shown in Fig. [5] 
leads to strong broadening and increase in height of the 
switching probability distribution seen in Fig. [2Ji. Be- 
cause motion near the saddle is diffusive, switching paths 
loose synchronization. In other words, the distribution 
of times spent by the system near the saddle point is 
comparatively broad. This is why it is advantageous to 
study the distribution of switching paths in the space of 
dynamical variables rather than in time. 



C. Lack of time reversal symmetry in a driven 
oscillator 

Another generic feature of the observed distribution is 
characteristic of systems far from thermal equilibrium. 
For equilibrium systems, the most probable fluctuational 
path uphill, i.e., from an attractor to the saddle point, 
is the time reversal of the fluctuation-free downhill path 
from the saddle point back to the attractor. More pre- 
cisely, it corresponds to the change of the sign of dissi- 
pation term in the equation of motion [45L loa |. i.e., to 
replacing T with — T in Eq. (|2ip . In ovcrdamped equi- 
librium systems with detailed balance, these two paths 
coincide in space (but are opposite in direction). 

Our parametric oscillator is driven far from thermal 
equilibrium. Therefore the uphill section of the MPSP 
does not simply relate to the deterministic trajectory 
with reversed sign of dissipation. This section of the 
MPSP, i.e., the most probable fluctuational path from 
A\ to the saddle point at the origin is plotted as the 
thick solid line in Fig. [7J Upon sign reversal of the dis- 
sipation, the attractor becomes a repeller, as in the case 
of systems in thermal equilibrium desribed earlier. How- 
ever, in contrast to equilibrium systems, it is also shifted 
away from its original location (from A\ to A[ in Fig. [7J. 
The dissipation-reversed path is shown as the thin solid 
line in Fig. [Jj In addition, Fig. [1] allows one to compare 
the uphill section of the MPSP with the deterministic 
downhill path from S to A±. Our data show that the 
uphill section of the MPSP, which is formed by fluctu- 
ations, the dissipation-reversed path, and the downhill 
noise-free path from the saddle to the stable state are all 
distinct. The time irreversibility of the switching paths 
is directly related to the lack of detailed balance of our 
driven oscillator, distinguishing it from bistable systems 
in thermal equilibrium. 




-12 -8 -4 
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FIG. 7: (Color online) Comparison of the MPSP and the 
dissipation-reversed path. The section of the most probable 
switching path from A\ to S is shown as a thick solid line. 
Upon changing the sign of dissipation, the attractor is shifted 
to a new location A[ and becomes a repeller. The fluctuation- 
free path with reversed dissipation from A[ to 5 is shown as 
the thin solid line. 
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VI. CONCLUSIONS 

In this paper we have studied the phase space distri- 
bution of paths followed in activated switching between 
coexisting stable states. The analysis refers to systems 
with several dynamical variables. We introduced a quan- 
tity, the switching probability distribution, that gives 
the probability density of passing a given point in phase 
space during switching. The distribution is defined in a 
way that makes it experimentally accessible. No a priori 
knowledge of the system dynamics is required except the 
positions of the stable states, which can be immediately 
determined, since the system spends near these states 
most of the time. 

The switching probability distribution was shown the- 
oretically to have a shape of a narrow ridge in phase 
space. Far from the stationary states, the cross-section 
of the ridge is Gaussian. The maximum of the ridge lies 
on the most probable switching path (MPSP). 

Experimental studies of the switching path distribu- 
tion were done using a high-Q micromechanical torsional 
oscillator. All parameters of the oscillator, including the 
nonlincarity constant, were directly measured. The oscil- 
lator was driven into parametric resonance, where it had 
two coexisting vibrational states that differ in phase. The 
paths followed in switching between these states were ac- 
cumulated and their distribution in the space of the two 
dynamical variables (the oscillation quadratures) was ob- 
tained. It was found that the distribution has indeed the 
shape of a Gaussian ridge. 

There is excellent agreement between the experimental 



and theoretical results, with no adjustable parameters. 
The measured maximum of the switching path distribu- 
tion lies on top of the theoretically calculated MPSP. 
The measured velocity of motion along the MPSP as a 
function of the position on the MPSP also quantitatively 
agrees with the theory. An important property of the 
path distribution is the total current conservation: the 
product of the velocity of motion along the MPSP and 
the cross-section area of the path distribution remains 
constant. This conservation of probability current was 
demonstrated experimentally. In addition, we observed, 
for the first time, that the lack of detailed balance leads 
to the difference between the uphill section of the MPSP 
and the noise-free path with reversed sign of dissipation. 

The observation of the most probable switching path 
reported here provides, in some respects, an experimen- 
tal basis for the broadly used concept of a reaction co- 
ordinate, which can be associated with the coordinate 
along this path. Our method does not rely on the specific 
model of the fluctuating system but only on the charac- 
teristics accessible to direct measurement. It applies to 
systems far from thermal equilibrium as well as to equilib- 
rium systems. Measuring the switching trajectories can 
help to determine the model globally, far from the stable 
states. It can also provide an efficient way of controlling 
the switching rates by affecting the system locally on the 
most probable switching path. 

This research was supported in part by NSF DMR- 
0645448 (HBC) and NSF PHY-0555346 and ARO 
W911NF-06-1-0324 (MID). 
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